#!/bin/csh 
setenv LC_ALL C
setenv LC_CTYPE C
setenv LC_COLLATE C
setenv LC_TIME C
setenv LC_NUMERIC C
setenv LC_MONETARY C
setenv LANG C

gmtdefaults -D > .gmtdefaults
gmtdefaults -D > .gmtdefaults4

#gmt set PS_MEDIA 19cx50c



#######################################################
set line = 25
set ps = 25_normalfault1.ps
set origin_x_min = 9000 #set min_cmp num
set origin_x_max = 12000 #set max_cmp num
set origin_d_min = 0.41 #set 
set origin_d_max = 10.5
set trench_cmp=7700
#set width = 38 # set figure width here
set vert_exa = 1 #set fig vertical exaggeration
set height = `echo "scale=2; ($origin_d_max - $origin_d_min)/3" | bc`
set cm = c
#######################################################

set dmig='../../segy/MC'$line'migkdinit/migkd_coffstk_wa_it0.sgy'


########################

set new_x_min = `echo "scale=2; $origin_x_min * 6.25 / 1000 " | bc`
set new_x_max = `echo "scale=2; $origin_x_max * 6.25 / 1000 " | bc`
set width = `echo "scale=2; $height * ($origin_x_max - $origin_x_min) * 6.25 *0.001  / (($origin_d_max - $origin_d_min)* $vert_exa)" | bc`
#set height = `echo "scale=2; $width * ($origin_d_max - $origin_d_min) * $vert_exa / (($origin_x_max - $origin_x_min) * 6.25 * 0.001)" | bc`
set dimension = `echo $width$cm/-$height$cm`
set fig_height = `echo "scale=2; $height+4" | bc`
set fig_width = `echo "scale=2; $width+4" | bc`
#set xylimits_km = `echo $new_x_min/$new_x_max/$origin_d_min/$origin_d_max/`
set xylimits = `echo $origin_x_min/$origin_x_max/$origin_d_min/$origin_d_max/`
set trench_km=`echo "scale=2; $trench_cmp * 6.25 / 1000 " | bc`
set left_trench=`echo "scale=2; $new_x_min - $trench_km " | bc`
set right_trench=`echo "scale=2; $new_x_max - $trench_km " | bc`
set xylimits_km=`echo $left_trench/$right_trench/$origin_d_min/$origin_d_max` #6900*6.25/1000

##########################
gmt set PS_MEDIA `echo $fig_height$cm`x`echo $fig_width$cm`

# convert sgy to grd file = 1
set convgrd = 1
set tscal=-0.75

set cmps=0.0008
set cmps=0.0008
set cmps=0.0008
#set dscal=-1000/6.25/1000/5*0.0008
set dscal=-0.064

if ( $convgrd == 0) then
	
segyread tape=$dmig endian=0 | surange > tmp.asc
set cminm = `awk '{ if ( $1 == "cdp" )print $2}' < tmp.asc`
set cmaxm = `awk '{if ( $1 == "cdp" ) print $3}' < tmp.asc`
set ntraces = `awk '{if (NR == 1 ) print $1}' < tmp.asc`
set dzm = `awk '{if ( $1 == "dt" ) print $2/1000000}' < tmp.asc`
set dmaxm = `awk '{if ( $1 == "ns" ) print ($2-1)*'$dzm'}' < tmp.asc`

echo sufile ntraces cmin cmax: $ntraces $cminm $cmaxm

segyread tape=$dmig endian=0 | sushw key=counit a=1 | suchw key1=cdp key2=cdp | segyhdrs format=5 | segywrite tape=tmp.sgy conv=0 endian=0
/usr/bin/segy2grd tmp.sgy -G$dmig.grd -I1/$dzm -R$cminm/$cmaxm/0/$dmaxm -V -Sc -M$ntraces                

		
segyread tape=$veld endian=0 | surange > tmp.asc
set cminvd = `awk '{ if ( $1 == "cdp" )print $2}' < tmp.asc`
set cmaxvd = `awk '{if ( $1 == "cdp" ) print $3}' < tmp.asc`
set ntraces = `awk '{if (NR == 1 ) print $1}' < tmp.asc`
set dzvd = `awk '{if ( $1 == "dt" ) print $2/1000000}' < tmp.asc`
set dmaxvd = `awk '{if ( $1 == "ns" ) print ($2-1)*'$dzvd'}' < tmp.asc`
		
echo sufile ntraces cmin cmax: $ntraces $cminvd $cmaxvd
 
segyread tape=$veld endian=0 | sushw key=counit a=1 | suchw key1=cdp key2=cdp | segyhdrs format=5 | segywrite tape=tmp.sgy conv=0 endian=0
/usr/bin/segy2grd tmp.sgy -G$veld.grd -I1/$dzm -R$cminm/$cmaxm/0/$dmaxm -V -Sc -M$ntraces                


set cminm = `grdinfo -C $dmig.grd | tail -1 | awk '{ print $2}'`		
set cmaxm = `grdinfo -C $dmig.grd | tail -1 | awk '{ print $3}'`
set zmax=40.0
exit		 
endif


#::::::::::::: DEPTH MIGERATIONS :::::::::::::::::::
gmt makecpt -Cwgb.cpt -T-6000/6000 > wgb_25.cpt
gmt makecpt -Crainbow.cpt -T1400/7000 > rbow.cpt
gmt grdmath $dmig.grd 0.016 MUL = tmp.grd
gmt grdimage tmp.grd -Cwgb_25.cpt -R$xylimits -JX$dimension -B+n -M -K > $ps
gmt psbasemap -R$xylimits_km -JX$dimension -Bxf5a5+l"Distance to trench (km)" -BS -O -K \
--MAP_TICK_LENGTH_PRIMARY=-0.02i --FONT_ANNOT_PRIMARY=4p --FONT_LABEL=4p --MAP_LABEL_OFFSET=1p --FONT_TITLE=4p --MAP_FRAME_PEN=0.5p >> $ps
gmt psbasemap -R$xylimits -JX$dimension -Bxf100a1000+l"CMP number" -BN -O -K \
--MAP_TICK_LENGTH_PRIMARY=-0.02i --FONT_ANNOT_PRIMARY=4p --FONT_LABEL=4p --MAP_LABEL_OFFSET=1p --FONT_TITLE=4p --MAP_FRAME_PEN=0.5p >> $ps
gmt psbasemap -R -J -Byf0.5a1.0+l"Depth (km)" -BWe -O -K \
--MAP_TICK_LENGTH_PRIMARY=-0.02i --FONT_ANNOT_PRIMARY=4p --FONT_LABEL=4p --MAP_LABEL_OFFSET=1p --FONT_TITLE=4p --MAP_FRAME_PEN=0.5p >> $ps
psbasemap -R -J -B+n -O >> $ps
psconvert -TG -P $ps -E1000
rm -f gmt.*
rm -f *grd
rm -f *sgy
rm -f *ps
